home *** CD-ROM | disk | FTP | other *** search
/ Amiga Collections: Camelot / Camelot 098 (1990-12)(Swedish User Group of Amiga)(SE)(PD)[WB].zip / Camelot 098 (1990-12)(Swedish User Group of Amiga)(SE)(PD)[WB].adf / XLisp-Stat / Functions / robustregression.lsp < prev    next >
Lisp/Scheme  |  1990-10-11  |  2KB  |  43 lines
  1. ; book pp.174-176
  2.  
  3. (defun robust-coefs (x y wf &key (weights (repeat 1 (length y)))
  4.                                  (tol .0001)
  5.                                  (count-limit 20))
  6.   (let ((x (if (matrixp x) x (apply #'bind-columns x))))
  7.     (labels ((as-list (x) (coerce (compound-data-seq x) 'list))
  8.              (rel-err (x y)
  9.                (mean (/ (abs (- x y)) (+ 1 (abs x)))))
  10.              (reg-coefs (weights)
  11.                (let* ((m (make-sweep-matrix x y weights))
  12.                       (p (array-dimension x 1)))
  13.                  (as-list
  14.                      (select (first (sweep-operator m (iseq 1 p)))
  15.                              (1+ p)
  16.                              (iseq 0 p)))))
  17.              (fitvals (beta)
  18.                (+ (first beta)(matmult x (rest beta))))
  19.              (improve-guess (beta)
  20.                (let* ((resids (- y (fitvals beta)))
  21.                       (scale (/ (median (abs resids)) .6745))
  22.                       (wts (funcall wf (/ resids scale))))
  23.                  (reg-coefs wts)))
  24.              (good-enough-p (last beta count)
  25.                (if (> count count-limit)
  26.                    (format t "Iteration limit exceeded ~%"))
  27.                (or (> count count-limit)
  28.                    (and last (< (rel-err beta last) tol)))))
  29.     (do ((last nil beta)
  30.          (count 0 (+ count 1))
  31.          (beta (reg-coefs weights)(improve-guess beta)))
  32.         ((good-enough-p last beta count) beta)))))
  33.  
  34. (defun make-wf (name &optional (k (case name (biweight 4.685)
  35.                                              (cauchy 2.385)
  36.                                              (huber 1.345))))
  37.   #'(lambda (r)
  38.     (let ((u (abs (/ r k))))
  39.       (case name
  40.             (biweight (^ (- 1 (^ (pmin u 1) 2)) 2))
  41.             (cauchy (/ 1 (+ 1 (^ u 2))))
  42.             (huber (/ 1 (pmax u 1)))))))
  43.